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Abstract 

In this paper we consider the numerical solution of the Hamiltonian wave equation in 
two spatial dimensions. We construct a two step procedure in which we first discretize the 
space by the Mimetic Finite Difference (MFD) method and then we employ a standard 
symplectic scheme to integrate the semi-discrete Hamiltonian system derived. The main 
characteristic of the MFD methods, when applied to stationary problems, is to mimic 
important properties of the continuous system. This approach yields a full numerical 
procedure suitable to integrate Hamiltonian problems. A complete theoretical analysis of 
the method and some numerical simulations are developed in the paper. 


1 Introduction 

Because of the symplectic structures, Hamiltonian partial differential equations (PDEs) are 
used to give a mathematical representation of many physical systems and are of interest to 
various applicative fields, see for instance quantum field theory, meteorology, nonlinear optics, 
weather forecast. 

An important requirement that any numerical method for Hamiltonian PDEs has to satisfy 
is the preservation of the intrinsic geometric properties of the original continuous problem. In 
particular, the numerical procedure should preserve the symplectic structure of the Hamiltonian 
system during numerical simulations. A standard procedure to derive a suitable method for 
an infinite-dimensional Hamiltonian PDE consists into two steps: in the first one the system 
is discretized in space in order to obtain a finite-dimensional Hamiltonian system, and then 
the semi-discretized system is solved in time by a symplectic integrator [21 Ha na mini HOI- 
There exists also a recent approach in which the space and time are considered on equal footing, 
this approach requires a multi-symplectic formulation of the system and leads to the multi- 
symplectic numerical schemes for the numerical solution of the PDEs (see miiadTiisD])- 

The effectiveness of this approach is ensured by the property that the derived semi-discrete 
system is a finite-dimensional Hamiltonian system of ordinary differential equations (ODEs). 
The space discretization of a Hamiltonian system is usually performed by one of the following 
techniques: finite difference methods, finite element methods, spectral methods, pseudospectral 
methods, Fourier expansion, wavelet based methods (see for instance (21 EZl HSl IH EH])- 
However, these semi-discretization approaches could become very expensive or could not be 
applicable when the space dimension d is greater than d = 1. 
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Instead, in this paper we consider the Mimetic Finite Difference (MFD) method to ap¬ 
proximate the continuous problem combined with a standard symplectic integration in time to 
integrate the derived semi-discrete Hamiltonian system. 

The main results about MFD methods, for stationary problems, can be found in the recent 
book [7] and papers EHiiia where, in particular, the theoretical framework of the mimetic 
spaces and the discretization of the operators are introduced. Significative applications of 
MFD methods may be found for instance in [T31 [IHl El 0131 HI] • Among the first publication 
in this field it is worth mentioning [311331 where a first approach to mimetic discretization of 
the continuous operators can be found and the fundamental papers mill] where the modern 
approach to MFD was introduced. A generalization of the MFD methods has been recently 
proposed, the virtual element methods (VEMs); we cite [11 (H |3T1 [3311331 [5] as a very short 
representative list. 

Recently in [23], MDF methods has been applied to the space discretization of PDFs of 
parabolic type in two dimension, showing how this technique preserves invariants of the solution 
better than classical space discretizations such as finite difference methods. 

The main characteristic of the MFD methods is to mimic important properties of the con¬ 
tinuous system, e.g., conservation laws, symmetry and positivity of the solutions, and the 
most important properties of the continuous differential operators, including duality and self¬ 
adjointness relations. Furthermore MFD methods can be applied for general polygonal and 
polyhedral meshes of the space domain instead of more standard triangular/quadrilateral grids. 

The main novelty of this paper is the use of MFD methods for the space discretization of 
the nonlinear wave equation in 2D coupled with a standard symplectic method (the implicit 
midpoint scheme) for the time integration. We derive a full numerical discretization procedure 
which will exploit the conservative properties of the MFD approach associated to the sym¬ 
plectic features of the time integrator. We show that the mimetic semi-discrete Hamiltonian 
is preserved in time and we derive the conservation law for the mimetic semi-discrete energy. 
Furthermore we give a bound for the conservation of the full discretized Hamiltonian and for the 
conservation of the full discretized energy. We also prove the convergence of the semi-discrete 
and fully discrete solutions to the solution of the original problem 

The paper is organized in the following way. In Section [2] we recall the basic elements of the 
MFD approach. In Section[3]we recall the mathematical form of the Hamiltonian PDF we wish 
to study. In Section 0] we apply the MFD method to the continuous problem and we give a 
result of the convergence of the semi-discrete solution to the continuous solution of the original 
problem; we define the semi-discrete Hamiltonian and energy density, show their conservation 
laws. In Section 13] we discretize the semi-discrete system by using a symplectic time integrator, 
the implicit midpoint rule, of the second order in time. We will prove the convergence of the 
full discrete numerical solution by providing an error estimate of the second order in space and 
time. Hence we give a result about the conservation of the discrete Hamiltonian and of the 
discrete energy of the system. Section [3] is devoted to show some numerical results. 

2 Background on Mimetic Finite Differences Methods 

In this section, for ease of reading, we recall the basic concepts and notations on MFD 
methods which will be used to discretize PDFs in the spatial domain H C where we assume 
D bounded polygon. For more details on this subject we refer the interested reader to the 
recent book [7] or to the papers [HHHlHn]. Let Lo a measurable subset of the domain H and let 
K e (L°“(H))^^^ a full symmetric positive definite tensor. By making use of standard notation, 
we consider the following scalar products: 



( 1 ) 


( 2 ) 


It is clear that, in the sense of distribution, 

(KVu, cr)]K_n = —(mj divcr)^ 2 (f 2 ) for all u G L^(H), cr £ iJ(div, H) 
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thus we get the duality relation with respect to the scalar product 0 and m 


KV = -(div)*. (3) 

Let Th be an unstructured mesh of into nonoverlapping simply-connected polygons with 
flat faces, where 

h := sup diameter(c). 
cGTh 

Let 8h be the set of edges of the polygons in Th- We use the following notations for the 
mesh objects: c & Th denotes a general cell in the mesh with measure |c| and centroid Xc; 
f € £h denotes a general edge of the cell c with measure |/| and centroid x/; n/ indicates the 
unit normal vector to the edge / with preassigned direction; acj = ±1 represents the mutual 
orientation of the vector ny and the outward normal vector to / with respect to the cell c. 

Moreover, let Th = Th,8h, and let a = c,f, then we denote with yh{<y) the subset of yh 
of all the elements that are related with cr, and we indicate with | 3 ^h(tT)| the cardinality of this 
set. For example Th{f) denotes all cells sharing face / and £h{c) denotes all faces forming the 
boundary of cell c. 

In the following we take on the element c G Th the shape regularity assumptions listed, for 
instance, in [ziiia. A possibility is to assume that for all h, each element c in Th satisfies: 

(Ml) c is star-shaped with respect to a ball of radius greater then 7 he, 

(M2) any two vertexes in c are at least cr hx apart, 

where he is the diameter of c. The constants 7 and a are positive and uniform with respect to 
the mesh family. 

The mesh objects will define the degrees of freedom of the discrete system, that is these will 
define the space of the discrete pressures and discrete fluxes. 

Let 

Ne:=\Th\, Nf:=\£h\, N* := ma,-x\£hic)\. 

C 

Let Ch be the set of the pressures that are piecewise constant on Th, i.e. 

Ch '■= {u G I u\e = const, Vc S Th } . 

Given a pressure u G L^(n), we define the interpolant discrete pressure G Ch with 

u dc, for all c S 7/i. 

The space Th of the discrete velocities is defined as follows. For all edge f G £h associate a 
real number Uf and we denote with cjh the vector with components given by the collection of 
all the {u}f}f^Sh- The symbol Th will represent the vector space of all ujh- Let a; G H{div,c) 
a vector function, and let us assume that all face-integrals 



J w • u/ dA, 


for all f G £h 


exist. Then the interpolant discrete flux of u in the space Th is defined by 
with 



w • n/ dS, 


for all f G £h- 


{^f)feSh 


Remark 2.1. The discrete spaces Ch, Th and the interpolation operators are defined starting 
from the degrees of freedom: 


• j^ude, for all c G Th, and u G L‘^{Cl), 

• w • n/ d^, for all / G £h(c), and cj G JI(div, c). 
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Remark 2.2. There are obvious correspondences: 


ur^{uc)cerh, and ui ^ {ujf)f^e^. 

With a slight abuse of notation we can refer to a function in the discrete functional spaces as 
a vector and vice versa. 

The definition of the mimetic scheme carries on with the discretisation of the differential 
operators. Let a; S H{div, c) with c gTh, then the Divergence Theorem states that 

/ diva; dx= a; • n dS' , 

J c J dc 

where n is the unit outward normal to dc. Therefore, the continuous operator div admits the 
immediate discretisation TJIV: Th —>■ Ch, with 

(T>XVa;/i)c = ^ acj\f\ojf 'i ujh & J^h- 

f&Shic) 

The operator VJV is called discrete primary operator. 

The next step in the construction of the MFD method is the definition of suitable inner prod¬ 
ucts on the discrete functional spaces Ch and Jtj that allow to construct the derived operators 
imposing the duality relations for the discrete operators. 

We assume, for the moment, the following scalar products on the vector spaces Ch and Th'- 

[uh, Vh]ch ■= u^Mc^Vh for all Uh,Vh eCh, (4) 

[u:h,eTh]j^^ ■■= uj^Mj^^cTh for all cjh, (Th & (5) 

where Me,, € e are suitable symmetric positive definite matrices. These 

matrices are locally constructed in such a way, on each cell, the corresponding local discrete 
inner products have to “mimic” the scalar products defined in ([T])) and ([2]). Therefore we would 
like that 


[uh,c,Vh,c]Ch,a =■ {uh,c)'^M.c^^^Vh,c ~ (uh, Vh)L'^(c), for all Uh,Vh e Ch, 

[i^h,c, crh,c]j^h,c =■ i^h,c)'^Mj^h,c(^h,c ~ (Wft , crh)K,c, for all UJh, (Th S Th, 

where, in general, with the notation rh^c we denote the vector with the degrees of freedom of 
the function r relative to the cell c. 

As regards the first local inner products, we observe that the vector rh^c has a single com¬ 
ponent, representing the (constant) value of Vh in the cell c. Then the only possible quadrature 
formula is 

['^hjC,'^h^c\Ch,c — — lolWcUc, 

therefore Mc^ ^ = |c| and Me,, := diag(|ci |,..., Icat,, |). It is clear that the discrete inner products 
gives the exact value of the continuous one whenever Uh,Vh S Ch- 

The definition of the local scalar product for the fluxes requires a different approach. The 
key idea is to define suitable consistency and stability constraints in order to introduce algebraic 
conditions on the elements of the matrix M_F;, Without spelling things out, we requires that 
the following properties are satisfied 

• consistency: let a;, cr two vector fields and let a;^, cr^ G J^h their interpolant functions. If 
(jj is constant in c and for each edge in f G £h{c), cr ■ is constant, then 


[<^h,c, crh,(^j^h,c 


K ^o; • (T dc; 


• stability: there exist two positive h-independent constants C* and C* such that 
C.,\c\{uh,c)'^i^h,c < {<^h,c)^J^h.c^h,c < C*\c\{u}l^^)uh,c £ J^h- 
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The last preliminary step in the construction of the MFD method is the dehnition of the 
derived discrete operators, which are obtained through a duality relation from the primary 
operators. Let us consider the spaces Ch, J^h equipped respectively with the scalar products 
From continuous duality relations m, we can introduce the discrete operator 

QTZAV: Ch ^ 


and impose the duality relation: 

[^h,GTi-A'Duh]j^^ = —[VIVu:h,Uh]ch (^hMj^i^GTZAVut = -Uh'DIV'^Mc^Uh , 

for all (jJh & ^h^Uh G Ch, from which it follows that 

gUAV ■= 


Finally we can introduce the discrete counterpart of the continuous operator div KV, by defining 
the operator 

'■ Ch ^ Ch 


given by 


Ah := VIVgUAV. 


( 6 ) 


3 The continuous problem 

Let C be a bounded polygon and let us consider the nonlinear wave equation with 
homogeneous boundary value problem 

Utt{x,t) = divIKVM(a;, t) — f'{u{x,t)) in D x (0,T) 
u(x, 0) = uo(x) , Ut(x,0) = vo(x) in D (7) 

u{x,t) = 0 on dCl X (0,T) 

where K G is a full symmetric positive definite tensor, and the source term /' is the 

derivative of a smooth function /: R R. We would observe that no particularly restrictive 
assumptions on /' are required, for instance f in the sine-Gordon equation or the ones of 
polynomial type with respect to u may be considered. For seek of simplicity we consider in the 
proof /' global Lipschitz, however the convergence results are still valid for /' local Lipschitz 
(see Remark 14.21) 

o admits the equivalent formulation 

ut{x,t) = v{x,t) 

Vt{xA) = divKVM(a;,t) — f'{u{xA)) 
where the initial and boundary conditions are given by 
u(x, 0) = uo(x), v(x,0)=vo(x) in D u(x,t)=0, 

I® is said Hamiltonian formulation of ([7]) for which the Hamiltonian 

'H.[u,v\:= J + ivu • KVu +/(m)^ da; (9) 

is invariant with respect to time t along the solution, that is 

j^n[u,v] = o. ( 10 ) 

The energy density of the system is defined by 

E{u,v):=^v^ + ^S/u-Kyu + f{u). ( 11 ) 


in O X (o,r) 
in O X (o,r) 

v(x,t) = 0 , on dfl X (0,T) . 
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The total derivative of E{u,v) with respect to t, along the solution {u,v) of ([HI), is given by 

Et = (divKVit)?; + Vu ■ KVu = div {v KVu). 

Let u}{u, v) := —vK.Vu the energy flux, then we have the energy conservation law 

Et{u,v) + divuj{u,v) = 0 , (12) 

which is more general than the global conservation of the Hamiltonian. Indeed if the energy 
conservation law holds, then it is easy to prove that -^'H[u,v] = 0. 


4 The semi-discrete problem 

By using the MFD approach we can approximate the continuous operators by discrete ones, 
in order to derive the semi-discrete problem for the wave O- Then the resulting semi-discrete 
wave equation reads: 


Uh,tt(t) = Ah Uh(t) - f{uh(t)) for t e (0, T), 

Uh(0) = Uh,o, Uh,t(0) = Vh,o , 


(13) 


where Uh,o ■= and Vh,o ■= Vq are the interpolant functions in Ch of the initial data. In the 
same way, m can be discretized in the following form 

( Uh,t{t) = Vh{t) forte(0,T), 

< Vh,t{t) = AhUh{t) - f'{uh{t)) fortG(0,T), (14) 

[ u/j( 0 ) = Vh{0) = Vh,o. 


We observe that the semi-discrete m preserves the Hamiltonian structure of ([S]). In light of 
the definition in Section [21 the Hamiltonian functional TL in admits the natural mimetic 
semi-discretization: 

'Hh[uh,Vh] ■■= ^[vh,vh]ch + ^[G'R.A'Duh,GTZAVuh]j^h + [fi'^h)A]ch : ( 15 ) 

that will be called mimetic semi-discrete Hamiltonian functional. 

We can observe now that, if we denote with the gradient with respect to the variable 
Vh and with the gradient with respect to the variable Uh, then 

Vv^'Hh[uh,Vh] = = Vh , 


and 


^cl'^ uh'^h[uh,Vh] = Mch {GE.AV^ 'Wlj^^G'R,AVuh + M.c^f'{uh)) = —AhUh + f'{uh). 

Hence (1131) may be written as a Hamiltonian system of ordinary differential equations (ODEs), 
that is as: 

where is the canonical symplectic matrix ( _ ) while V denotes the gradient with 

respect the variables Uh and Vh- 

We can conclude that the MFD approach gives a finite-dimensional system of ODEs that 
retains the Hamiltonian character of the given PDE. Therefore MFD methods can be considered 
powerful scheme for the spatial discretization of Hamiltonian PDEs. 
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4.1 Convergence for the semi-discrete problem 

Now we will investigate the convergence of the solution Uh of the semi-discrete wave equation 
m to the solution m of ([5|) in norm. Before analysing the error between the solution, 

we have to show some preliminary technical results. 

Let us introduce the energy projection Vh ■ —>■ Ch, with u ^ Vtu defined as the 

solution of the diffusion problem 

{ find VhU € Ch and cr^ £ Th such that 
[(Th, ^h]Fh + [Ph U, VXV 0Jh]Ch = 0 for all G (16) 

[VIV (Th, Wh]ch = [(divKVu)^, Wh]ch for all Wh G Ch- 

In particular, for the duality relation between the operators, the projection Vh u satisfies 

AhiVhu) = {divKS/uY. (17) 

In the following we use || • \\c^ to denote the norm induced by the scalar product [•, 
(that is equivalent to L^(n) norm on Ch)- We will denote with C a generic constant, possibly 
different at each occurrence, independent from the mesh size h and the time step size t. In 
order to prove the convergence results, we need the following Lemma (see m for the proof). 

Lemma 4.1. Let us assume the convexity of the domain Cl. Let u £ and let VhU the 

energy projection of u. Then the following estimate holds: 

\\u^ - Phu\\c^ < Ch'^\u\H2(n). 

While the next Lemma shows the spectral properties of the operator A;i(see [5S] for more 
details). 

Lemma 4.2. The spectrum a{—Ah) of —Ah satisfies 

(T{-Ah) C [s„ s*h-% (18) 

where s» and s* are positive and h-independent constants. 

For the treatment of the nonlinear term we have the following lemma. 

Lemma 4.3. Let f be a smooth function and let u £ Then 

WfiuY -nuY\\c,<ciu)h^ 

Proof Let Uc '■= {u^)c for every c € 7h- Then the nonlinear term may be treated in the 
following way. Using the Taylor expansion, since u £ H^{Lt) 

f'{u{x)) = f'{uc) + f"{uc){u{x) - Uc) + ^f"'(Ucix)){u{x) - Mc)^ for a.e. a; £ c (19) 

for suitable Uc{x) and for every c £ Th- Then, setting f{u)c '■= [f{u)^)^ and using m , we 
have: 

f{u)c - f'{Uc) = n / {f'{u{x)) - f'{Uc)) dx 

^ R / f"{uc){u{x) - Uc) + {ucix)){u{x) - Uc)"^ - /'(mc)^ dx. 

Now, since f"{uc) is constant and, by definition, Mc = R fc'^(^) obtain 

f'(u)c - f'{uc) = (t{u)c 

with 

:= [ f"'(uc{x)){u{x)-Uc)^ dx. (20) 

2 c Jc 


7 


Now we observe that u G H^{n) for classic Sobolev embedding theory, implies that u G 
then, since Uc{x) is bounded by u(x) and the constant value Uc we obtain that Uc G L°°{c) for 
all element c G Th- Therefore being / a smooth function, f"'{uc) G L°°{c). Now, using the 
Holder Theorem 


{u,t)\\c^ f f"iucix)){u{x) -Ucfdx <\Y i l/'"(«c(a;))l(M(a:) - u^f dx 

cdTh cdTh 


< ^ E iir(«c)iUoo(c)ii(«-«c)^iiLi(c) <c E ii(u-uc)"iui(c). 


ceTh 


c&Th 


Now, using standard polynomial approximation results [U], we have 

\\cr{u,t)\\c^ <cY II(«-^^c)^||l 1 (c) = C* E ll(^-'“c)lli=(c) 

ceTh c^Th 

< C E C 

cGTh 


( 21 ) 


( 22 ) 


□ 


Now we have the instruments for proving the following convergence theorem. 

Theorem 4.1. Under the assumptions of Lemma \4-1\ and Lemma \4-‘A u{x, t) he the solution 
of m and Uh(t) be the solution of (fTl?ll . Let us assume that u{-,t) G for all t G [0,T]. 

Moreover let us assume that f is globally Lipschitz. Then, for all t G [0,T], it follows that: 


\\u{t)^ — Uh{t)\\ch < Ch? (|uo|i/ 2 (Q) + |uo|i/ 2 (n) + \ut\L i(o,t,ff=(n)) + 

+ \utt{t)\L^(o,t.m{Q.)) + \uii)\L^(o,t,min)) + \u{t)\h{o,t,m{n))) > 


where u{tY denotes the interpolant ofu{x,t) in Ch and the scalar function ifif) is bounded for 
all t G [0, T], 

Proof. The proof follows the guidelines of Theorem 1 in [5] for given for the finite element 
approximation. Let us set 

Uhit) - u{tY = {uh{t) - Vh u{t)) + {Vh u{t) - u{tY) =■■ 'd(t) + g{t). (23) 


We study separately the two terms. The second term represents the error generated by the 
energy projection; using Lemma l4.11 we obtain 

ll£'(0llc;, = W'PhuY) -uYYWc,^ < Ch‘^\u{t)\H2(a) =Ch^ (^\u{0)\H'^(n) + J \ut{s)\H2(n)ds'^ 

< Ch^ (|uo|H2(n) + |Mt|Li(o,t,//2(n))) • 

(24) 


For the first term, from (HSl) and CH), we get 

duit) - Ah 'd{t) = -f'{uh{t)) - {Vh u{t))tt + {diYKyu{t)Y 
and, since u is the solution of m, we obtain 

-dtt{t) - Ah -dit) = -f'{uh{t)) - Vh uu{t) + UuYY + {f'{u{t))Y 
= -Quit) - {f{uh{t)) - f{u{t)Y) 


and in particular 

[dtt{t), xIch - [Aft d{t), xW = -Ytt{t), x\ch - [f'{uh{t)) - f'{u{t)Y, x]ch (25) 




for all X S C/j. For t G [0, T] let use define 


G{t) := f -(/'K(s)) - fHs)y) ds, (26) 

Jo 

and let x = x(^) ^ in ((^51) be a function of t. Then it is straightforward to see that 

- Xt{t)]cH - [Ahd(<), x(t)W = - Uh,t){t) + G{t), xit)]cH + 

+ Xt{t)]ch - {G{t), Xt{t)]ch- ( 27 ) 


Let us fix r G [0, T) and we set in ((77)1 

X{t) ■= ^(s) ds, for t G [0, T], 

in particular we can observe that Xt{t) = ~^(i)- Now the duality relation among discrete 
operators and simple computations yield 


[Mt), '&it)]cH - [ST^-A'^Xt{t), GTZAVxit)]j-H = ^[(4 - Uh,t){t) + G{t), xit)W + 

- d{t)]c„ + [G{t), d(t)]c^ 

thus 

-UHM) + G{t), x{t)]c, + [G{t) - mW- 

(28) 

Integrating (EHl) with respect to t from 0 to r, observing that x(t) = 0, G(0) = 0 and by 
definition Uh^t{0) = ut{0y, we get 

ll^(r)||^, - 11^(0)11^, + ||I/7^^77x(0)||?., = 2 r[G(t) - gyt), mW dt 

Jo 

and then 

mr)\\l < 11^(0)11^, +2 r||G(t)||cJ|^(t)||c.dt + 2 T ||^.(t)||cj|^(t)||c.dt. (29) 

Jo Jo 

Now by definition (j26p . from Lipschitz assumption on the load f ^ and Lemma 14.31 we have 

||G(t)|k < r||/'(u,(s))-ms))^|kds 
^0 

< t\\r{uys))-r{u[sy)\\c,ds+ \\\r{u{s)y-r{u{sy)\\c,ds 
Jo Jo 

<G f \\'&is)\\c^ds + G f Wgisyc^ds + Gh"^ f |w(s)|^i(o) ds. 

Jo Jo Jo 

Therefore, from Cauchy-Swartz inequality and 

2 ri|G(t)||cJ|d(t)||c.dt 
^0 

<2G f (f \\iJis)\\c^ds+ f \\gis)\\c^ds + h^ f |u(s)|^i(f 2 ) dsh|^?(t)||c^ dt 
Jo \Jo Jo Jo / 

~^Io (/ dt + G {1 + T)\\^{t)\\l^ dt 

<G{u)Th^ +G{1 + T) rWmWk^t 
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and always from Cauchy-Swartz and (IMll 


2 r\\gmcjm\\c,dt< rMt)\\kdt+ r\m\\kdt<c{u)Th^+ rwmwkdt. 

Jo Jo Jo Jo 

By collecting the previous estimates in ([2^ . from (pi|) we get 

mT)\\k<\\m\\k + c{u)Th^ + {i+c+cT) rwmwkdt- m 

Jo 

It is straightforward to check that 

ll^(0)lk = ||^^(0)|k<C(uo)/l^ 


then by Gronwall inequality it holds that 

MT)\\k<Ciu,uo)h^Te^^. 

from which follows the thesis. 


□ 


Remark 4.1. The use of the projection Vh u in the proof of the theorem seems to be necessary. 
Indeed if we compute directly Uh{t) — u{t)^ as done for example in [3S], we obtain a term of 
the form 

L{u) := ||(divKVit)‘^ — AhU^k i 

II II C/i 

and L{u) does not converge to zero. For instance in Figure[T]we plot the asymptotic behaviour of 
L{u) as a function of h for u{x, y) = sin(7ra;) sin(7r?/), tensor K = /2 and domain Vl = [0,1] x [0,1] 
discretized with the sequence of Voronoi meshes introduced in Section^ see Figured) The value 
of L{u) does not seem to converge to zero as h is reduced. 



Figure 1: Asymptotic behaviour of L{u) as a function of h. 

Remark 4.2. Using standard theory of polynomial approximations and the definition of the 
projection Ph (see 13), it is possible to extend the proof of Theorem 14.11 to the case of /' local 
Lipschitz instead of global Lipschitz. The proof being analogous but more technical. 

4.2 Conservation laws for the semi-discrete problem 

As for the continuous system, it is easy to prove that the global semi-discrete conser¬ 
vation law of the Hamiltonian semi-discrete functional RhluhjVh] is preserved. Indeed using 
the duality definitions of the discrete operators, we have 

■^Phluh,Vh] = ^ {^[vh,vh]ch + Uh.QPAVUh]Th + [fiuh)A]ch'^ 

= [vh,t,vh]ch + [GP-AV Uh,t,GTlA'D Uh]j^h + [f'iuh)uh,tA]ch ( 31 ) 

= [AhUh,Vh]ch - [f'{uh),Vh]ch + [GPAVvh,GPAVuh]j^,, + [f{uh),Vh]ch 
= —[GPAVuh,GPAVvh\j^f^ + [GPAVVfnGPAVUh\j^^ = 0 , 
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along the solution {uh(t),Vh[t)) of (fTlTll . 

We can define the mimetic semi-discrete energy density Eh S Ch with 

Eh{uh,Vh)\^ ■.= ]^\c\vh^^A-]^[{QTlA'D Uh)cAQ'E,AV Uh)c]jr^^^ + \c\ f{uh,c) , ( 32 ) 

and by computing its derivative with respect to t along the solution, we have 


^Eh{uh,Vh)\, = ^ {^\c\vl^c + \ [{QTIAV Uh)c,{Q'R,AV uh)c]j^^_^ + |c| 


= |c| Vh,c -T:'Vh,c + 
at 


dt 


{QTZAV Uh)cj {GTIAV Uh)c 


Tk, 


+ |c| f'{uh,c)^Uh,c 
dt 


= |c|(A,j Uh)cVh,c - |c| f'{Uh,c)vh,c + [{G'R.A'Dvh)c, {G'R-A'Duh)c]jr^^ + |c| f'{uh,c)vh,c 
= |c| {Ah Uh)c Vh,c + [{GTIAV Vh)c, {GEAV Uh)c]jr^^ ■ 

Then, the following mimetic semi-discrete energy conservation law holds: 


Eh,t{'^h,Vh) + Eh{uh,Vh) = 0 (33) 

where Fh{uh,Vh) G Ch, defined by 

Eh,c{'^h,Vh) = -\c\{AhUh)cVh,c- [{GTlAVvh)c,{G'RAVuh)c]j7^^ , 

is a natural discretization of DIV {vhG'EA'Duh)- (13dll represents the mimetic approximation 
of the energy conservation law CH)- 

We have observed that the continuous Hamiltonian and semi-discrete Hamiltonian are first 
integrals respectively for system (|5]) and m- In particular, for all t G [0, T] we have 

'H[u{t), u(t)] = Hluo, uo] =: Ho, and Hh[uhit), Vh{t)] = H[uhfi, =: Hh,o 

where {u{t), v{t)) is the solution of ([5]) and {uh{t), Vh{t)) is the solution of (fTil) . In the following 
results we estimate the error between the continuous Hamiltonian and the semi-discretized 
Hamiltonian. 

Lemma 4.4. Let u G ffQ(H), and v G i?g (fl) and let and G Ch their respective interpolant 
function in Ch- Then it follows that: 

\H[u,v]-nh[u\v^]\ < Ch"^ (|w|Hi(n) + ^1^1(0) + |u|ffi(f2)|w|ff3(o)) (34) 


Proof. We split the bound for the three terms composing the Hamiltonian functional. Let us 
start with 


f{u{x))dx- [f{u^),l]c^ 


and we observe that since u G L^{Gt), using the same computations in (fTiHl cell by cell 


( 35 ) 


fiu{x)) = f{Uc) - f'{Uc){u{x) - Uc) + -f"{Uc{x)){u{x) - Ucf 


where we observe that using the same arguments in the proof of Lemma 021 the term f"{uc) G 
L°°{c). Now, since by definition Uc = -^ f^u{x)dx, using the same computation in (1^ and 
(03), it follows that 


/(u(a;))da: - [/(uj, l]c;,., 


(/(u(x)) - f{uc))dx 

= iucix)){u{x) - Ucf dx 


<Chf | m |^ i ( c )- 


11 












By adding in the cell c G Th bound the first term as follows 


f{u{x))dx- [f{u^)A]cH 




ceTh 


For the term 


/ v'^ dx - [v^,v^]c^ 
Jn 


( 36 ) 


we observe that 1 ]c^j thus we can use the computations used before 

with /(s) = s^. 

For the last term, we preliminary observe that, using integration by parts and homogeneous 
boundary conditions, and since QTZAT) = —T>IV*, we have to estimate 

/ (divIKVM)uda; — 

Jn 

Now, by definition of interpolation operator in Ch and (1171) . we get 

/ (divKVu) u da; — [A;, , u^]ch = 

Jn 

= j ^divKVu — (divKVu)^^ {u — u^) dx + J (divKVu) da; — [A/i , u^]ch 


JCh 


=: a + (divKVu)^ — A^i M'^, 

= a + [A,, {VhU- u^),u^ - Vhu]^^ + [A/i {VhU- u^),T’hu]^^ 
= a + [A/i {Vh u - u^),u^ -'Phu]^^+ VhU - , (divKVu)^ 

=: q; + /3 + 7 . 

Using standard polynomial approximation estimates [S], we have 


Ch 


a = 


< C |M|Hi(n) |w|_f/3(n)- 


J ^divKVu — (divKVu)^^ {u — u^) da; 

Moreover, by Lemma HIT] and Lemma we get 

1 ^ 1 = [/^h{VhU-u^),u^ -Vhu]^^ <\\Ah\\\W -Vhu\\c^ 

<Ch ^ |u|^ 2 (f 2 ) = C|u|^ 2 (Qp 

For the last term 7 , we it holds 


I 7 I = 


VhU — u^, (divKVu)'^ 


^Ch 


Finally, by collecting (1551) . ((5^ and (liUll in (1571) . we obtain 


(divKVu) u da; - [A/j Uh, Uh]ch 


Finally, the thesis follows from ((55|) . (I5S1) and (Hill . 


(37) 


(38) 


(39) 


<C\u\H^n) \\VhU-u^\\c^ < Ch^ |u|ff2(n)- (40) 


<ChJ + |u|//i(o) |u|//3(Q)^ . (41) 


□ 


Since in order to define the initial data in the MFD discretization we used the interpolantion 
operator in Ch, Lemma 14.41 implies the following estimates on the error between the continuous 
and semi-discrete Flamiltonian. 

Theorem 4.2. Let (u(t),v(t)) be the solution of system ([§]) and let (uh{t),Vh(t)) be the solution 
of (1141) . Then, for all t G [0, T] it holds 


\'H[u{t),v{t)] - 'Hh[uhit),Vhit)] \ < 

where the constant C depends only on the regularity of the initial data uq and vq ■ 


(42) 
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Proof. For all t G [0;T], it holds that 

'H[u{t),v(t)] =no, and 'Hh[uhit),Vhit)] = V-hfi- 
Now, since Uh,o = Wq and Vh,o = Vq, from Lemma [4.41 we get the thesis. □ 


5 The fully discrete problem: a symplectic MFD method 


In this section we will derive a symplectic mimetic finite difference method by ap¬ 
plying a time integrator scheme to the semi-discrete problem (US. Because of the Hamiltonian 
structure of the system (1131) a symplectic scheme is usually employed to integrate in time, in 
order to preserve the symplectic structure of the flow map of the system. 

Thus, we apply the symplectic implicit midpoint (SIM) (which is a scheme of second 
order in time) to problem p3p and get: 


,,n+l 


7/”+l — 7/" 4- T-- 
U), — Uh + T- 




. '^h — Vh = '^h,0 


,,n+l 


-f 


,,n+l 


(43) 


or equivalently 


n+1 


.."U+l 


.."U+l 


-f'[ul+T- 


n+1 


— Uu + T- 




= Uhfi, Vf^ = Vh,o , 


where denotes the numerical approximation of {uh{t), Vh{t)) at time = nr, for 

n = 0,...,N and r = T/N represents the time step length. Finally, by eliminating and 
, we can express the system in the following form: 




,,n+l 


2t 

= 2 (A. 


,U+1 


2t 

2ui: 


(/'(«r*)+/'(»:")) 


(44) 


where we use the notation 


.9+5 


,9+1 


+ < 


ul " = 


for q = 0,..., TV — 1. 


5.1 Convergence for the fnlly discrete problem 

We investigate the convergence of the sequence {u'f}n=i,....N to the exact solution u of 
problem ©• The following result states the convergence of the numerical procedure in discrete 
norm. 

Theorem 5.1. Let u he the solution of problem © and let he the sequence 

generated by Then, if u G C^{[0,T], H'^{Ll)), it follows that: 

W^itnY — w)(||ch < C h^ (|uo|ff2(n) + 9/2(0))) -I- C(45) 

where the constant C depends on the regularity of u. 
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Proof. Let us split the discrete error in the usual form 


Ul - uitnY = « - Vh U{tn)) + {Vh U{tn) - uitnY) = O’” + g”. (46) 

From Lemma l4.11 using the same argument in (1241) . we get 

lle^llc;. = \\rhu{tn)-u{tnY\\c^ < C'(|uo|(n) + jut|fl-i(o,t„,_f/ 2 (n))). (47) 

The analysis of the term cr" is more involved. We start by considering the first time step n = 0 
and we observe that using again Lemma 14.11 it holds that 

\W°\\Ch = llWft - 'Ph u(0)||Ch = \\uh,0 - Ph UollCh < |Mo|i/2(Q). (48) 

Let us analyse the first time step ti = t. Using the regularity assumptions on the solution u in 
the time variable, we have that for all a; G U it holds 

T 

u{x,t) = Uo + Tt;o(a;) + T t) + 

2 ^ (49) 

= Uo + TVoix) + Y (^divKV (^u [x, - /' [x, ^)) + -^ 

where R = O(t^) is the rest in the Taylor expansion of u{x, •). By definition (1481) with n = 0, 

ul = Uh,0 + T Vh,0 + Y {^h uY^ - f (uY^^ ) ■ 

Then, using CZD, recalling that u is the solution of © and interpolating (14911 in Ch it follows 
that: 

2 2 

2^1/2 - Y A;, ul/2 ^ ^ ^0 _ u{T/2Yj 

2 2 

= Uh,Q +TVh ,0 - “ Phu{T) + Y(divIKV {u (t/2)))^ + (7° 

2 2 

= Uh,0 + TVhfi - yY ~ Phu{T) + y(^« Y/Y + f (u (t/ 2)))^ + Cr° 

2 

= Uh ,0 + TVh ,0 + Y {Y - Y (YY^y) ~ Phuij) + (m(t) -uo-tvo- rY + cr° 

= Y {,Y {u{t/ 2)Y - f (uY^)) - (PhU(T)-u(TY)-R^ +(T°. 


Let us compute the scalar product of both sides of the previous equation with obtaining 


ai/2 ,1/2 


T 

Ch 2 L' 
2 


Ch 


(/' (u (t/2)Y - f - {Ph u{t) - uijY) -R^ + cr°, ( 7^/2 

1/2 


Ch 



u{t/2) -uY"^ 

+ 

Ch 

/^2 


^2 



a'/2 

T 



Ch 


\Ch 


^Ik+r^ 


Ch 

^ 1/2 


Ch 


(50) 


Now, since from Lemma SJ it follows — [AhVh,Vh]ch ^ 0 for all Vh, for small values of r using 
(|T7)1 and (gSl), we get 


ri/2 


Ch 


< C |Mo|ff 2 (n) + \ut\H^(o,T,H^(ci)) 
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and we can conclude that 


o’^llch < 2 ^ + ||o'°||cfe < ^|uo|ff2(n) + . 


(51) 


Now, we bound the error for a general time step n > 1. It is easy to see that the following 
relations hold 


l{tn+l) - 2u{tn) + U(tn-l) 


— “t" R , 


u{tn+l) + 2u{tn) + U{tn-l) = 4M(t„) + R , 

f'{u{tn+l)) + f'{u{tn)) 


f (u (in+l)) +R 


where R = O(r^) denotes the general rests in the Taylor expansion. Using the previous Taylor 
expansions, the definition of the scheme and (113, and recalling that u is the solution of 
©, we have 


cr"+l - 2cr” + cr”-l 1 . / + 2cr” + cr""! 

-5- 


(/' (C")+/' 


U{tn+l) - 2u{tn) + U{tn-l) 


{ divKV 


u{tn+l) + 2u{tn) + U{tn-l) \ V 


4(/'(«r*)+/'(«r’)) 


/ “t” 2utt{i’n) T — 

V 4 


+ 


f{u{t„+i)) + 2f'{u{tn)) + f'{u{tn-l)) 


+ (52) 


-Vh 


u(tn+l) - 2u{tn) + M(tn-l) 


= i (/' {u - /' (C^)) + ^ (/' (« {t. 

+ (uuitnY — 'PhUttitn)) + R 
= a'^+i +a^-i -g^^ + R^ 




+ 


where a? = 4'(“(U))^ f'iO ^ ^ ± i ^ 

Now, let 


_n+l _n 

5"+^ := ^ 


and let us observe that the following relations hold: 
cr"+l _ 2cr" + cr"-l ^"+i _ 


,^1 1 cr”+5-cr”"2 

S^+2 + = 2 - 


(53) 

I ~ I I 

Let us make the inner product of both sides of (I52|l with S^~^i + S'^~Y For the first term of 
the left-hand side, using (1^ . we get 




JCfe T 




Ch 




Ch 


(54) 


For the second term of the left-hand side in (EH), using (1^ and since A/j is self-adjoint, we 
have 




-r'^+h 


Ch 


Ch 


hcr ^ ,a 


Ch 


(55) 
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To bound the right-hand side, we preliminary observe that using the same argument of the 
proof of Theorem 14.11 and Lemma 14.41 it holds that 

ii«’iic. = ^ - /'K)|ic^ < c Unity - uiw,^ < c (ik^ii,^ + . 

Therefore, using the previous bound, the Cauchy-Schwartz inequality and the usual estimate 
in p, we derive 


:"+2 + a ^-2 + 5’^-2 


JCfe 


< C 

< c 






Ch 

2 

Ch 


Ch 

2 

Ch 


Ch 

2 

Ch 


+ T‘ 


S^+i +5^-^ 


Ch 


S^+i 


Ch 


+ T- 


<c( 

cr"+5 

2 

+ 

^'1^— ^ 
G 2 

2 

-b 


2 

+ 


V 


Ch 


Ch 


Ch 



Ch 


+ llt’ttllc^ + 


S^+2 +S^-2 


Ch 


+ T 


< C CT 


n+1 


ic.+2iKii?,+iK-‘ii;. 


j>‘+4 


Ch 


S^-2 


Ch 


-T^ + h^ 


Collecting (ISil) . (1^ and (1^ in (1^ . we obtain 


(56) 


- ( 

S^+i 

2 

S^-i 

2 

Ah 

-b 

Aha^ 

T \ 


Ch 


Ch 


Ch 



Ch 


< c \w 


n+1 


Ic,+2|k lie,+ lr lie. 


Moreover, (1^ and some simple calculations give: 

1 
T 


r+5 


Ch 


+ 


5^-2 


Ch 


-f + -b + . (57) 


IlCfcl 


- lk”“^llcj = 


<C[\\a-^yi,+2\\ay\l^ + \\a--yi^ 




Ch 


S^-2 


Ch 


(58) 


Now let us define 

pn __ 


S^+2 


Ch 


k”+'llc, + lk""" 


Ch 


Ah ct"+=,ct”+2 


Ch 


Using the estimates (1371) and (1351) . recalling that the operator — Ah is positive definite, we derive 
that 

pn pn—1 

- < C (+ -b r2)2 + C (T’" -b r”-i) , 

T 

and, by using the discrete Gronwall inequality, we obtain 

T’" < ^r° -b ^ T (/i^ -b r2)2^ . 

Now, using analogous arguments in ((3D1) and recalling bounds (1351) . ((3T]) . we obtain 

r° < C'(+ -br2)2, 


and thus 

lk”|ll <r'^<C(/i2+r2)2e"4Ct„, 

Hence, since tn < T, the above bounds gives 

\Wy\ch<Cih^+ry , (59) 

for all n = 1,..., A^, and collecting (llTll and (1331) in (1331) we get the thesis. □ 
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5.2 Conservation laws for the fully discrete problem 


The following result shows how the fully discrete method, built combining the MFD method 
and the symplectic implicit midpoint scheme, preserves, within an order of approximation, 
the Hamiltonian functional. Using classical results on the symplectic integrator methosd (see 
for instance m) we can state the following theorem about the long time stability of the Hamil¬ 
tonian. 


Theorem 5.2. Let be the sequence generated by system (H5|) . Then ifT < for 

a suitable positive constant 7 , it holds that: 


[Hh 






-Uh 


,0 .,01 


< Ct^ 


(60) 


Remark 5.1. It is well known that the SIM preserves the quadratic first integrals. Therefore 
if the load term / is quadratic, i.e. /(s) = fcs^, with k constant the Hamiltonian is exactly 
preserved along the solutions. 

By collecting the estimates of Theorem l4.2l and Theorem l5.2l we can provide a bound for the 
error in the Hamiltonian of the fully discrete procedure, stemming from the MFD discretization 
in space and the SIM integration in time. 

Theorem 5.3. Let {u{t),v{t)) be the solution of problem ([S|) and let be the sequence 

generated by system (1431) . Then, it follows that 

\nh[u^,v^] -'H[u{tN),v{tN)]\ <C{t'^ +h'^). (61) 

In Section [3] we have introduced the semi-discrete Energy density conservation law. Now we 
analyse the effect of time discretization in the semi-discrete Energy density conservation law. 

Theorem 5.4. Let be the sequence generated by system (SSI), and let for all c G Th 

and for all n 

EuAul.vl) := i|c| AAf + i [{gUAVulA + |c| /«J 


and 

Fn,c {uAvA = -|c| (A,<), vA - [(gnAVv ^)^, {gUAVulU^^^ . 
Then, the following estimate holds for all n 


Eh,cAA\v'^+^ 


h > ) ^h,c(u1,vff) U,"+5 

^ -I- rh,c 


< Clcl A . 


(62) 


Proof. We observe that, using SSI, it follows that 


1 


[[{gnAVuA^, {gnAVuA%]^^_^ - [{gnAVuiA iGnAVui),A^^) 

= (gizAVvl^^^ ,(gnAVul^^ 


^h,c 


and 



n+ 


1 V 

2 


n+1 

h,c 


— V 


n 

h,c 


T 


\,c 





Therefore, by collecting (1331) and (IMl) . we get 


(63) 


(64) 


EuAuT\vT") - EhAuAvA 


F, 


- G 


n+T 


n+7 


7Kr)-/Kc) 


n+' 

'A.c 


7'(C^)). (65) 
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Now, by the Taylor expansion, we derive 


/(C) - /«.c) = {utc - <c)f {ntj) + 

where R = 0{t^) denotes the rest, and thus 

+ Fh^c ( 


Rt'^ = T 


<:v' (C’) 


R 


EhAK^\K^^) - EhAK.vJ:) 


ra+= ra+'^ 


i) =i?|c 


□ 


6 Numerical tests 


In the present section we present some numerical results for the fully discrete case, i.e. SIM 
coupled with the MFD spatial discretization. The convergence of MFD has been evaluated in 
the discrete relative L^(n) norm of the difference between the interpolant A € Ch of the exact 
solution u and the numerical solution Uh at the final time T, i.e. 

p ^ \\A(T) - Uh,N\\cH 
■ hAAlA ■ 

Moreover we tested the total error in the Hamiltonian functional at the final step N, among 
the discrete solution and the continuous solution, that is: 

(Th,r ■■= lUhluAvh] - 'H[uo,vo]\ . 

We tested also the conservation of the Hamiltonian functional with respect to time integration, 
that is 

Sh,T ■■= \'Hh[uh,Vh] 

and the error in Energy density conservation law that is: 


Eh T '■= max 
’ c^Th 


Er,A<X)-EHAu^H-\vH-") 


T 


+ Fh^c 




We have considered the spatial domain ft = [0,1] x [0,1] C K^, and a general sequence of 
Voronoi meshes with h = 0.2,0.1,0.05,0.025 (see Figure!^, and r = 0.1,0.05,0.025,0.0125. 
For the generation of the Voronoi meshes we used the code Polymesher in . 



Figure 2: Sequence of Voronoi mesh with 0.05. 

Test 6.1. We consider problem ([5]) with the material tensor K = /2 and the load term / = 
A, where the initial data uq and wq are chosen in accordance with the exact solution 

u{t, xi,X 2 ) = sin(t) sin(7ra:i) sin(7ra;2)- (66) 

We implement the fully discrete problem in the time interval [0,1] with the SIM coupled with 
the MFD discretization for the sequence of polygonal meshes introduced above. In Table ?? 
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Table 1: Eh^r^ 5h,T and eh,T for fixed time step size r = 0.001. 



h = 0.2 

h = 0.1 

h = 0.05 

h = 0.025 

Eh,T 

7.2713323e-01 
7.5726618e - 03 
2.6233230e - 01 

1.9846010e-01 

1.9485290e-03 

1.4264664e-02 

5.2502301e-02 
5.0100939e - 04 
1.5776230e-03 

1.3086316e-02 

1.2559774e-04 

4.1846647e-04 




Figure 3: Asymptotic behaviour of Eh^r and 5h,T as a function of h for t = 0.001. 


we choose r = 0.001 and we show the errors in the solution Eh^r, for the Hamiltonian ah,T and 
for the Energy eh,T for different values of the mesh size h. 

In Figure [3] we plot the asymptotic behaviour of the errors in the solution and Hamiltonian 
as a function of h, in accordance with the theoretical order of convergence h'^ . 

In FigureHlwe plot the asymptotic behaviour of the errors in the Energy density conservation 
law eh,T at the final step N as a function of h ior t = h. We observe that, using (|ti2l) . we expect 
an order of convergence. 



Figure 4: Asymptotic behaviour of eh,T as a function of h for t = h. 

In Figure Owe extend the time interval setting T = 100 and we show the behaviour of the 
error 5h,T in the discrete Hamiltonian functional along the sequence with respect to 

the initial value for h = 0.05 and r = 0.001. We can observe that the Hamiltonian 

is numerically preserved by SIM. This results is in accordance with Remark 15.11 indeed in the 
test we are considering a quadratic function f{u). 

In Figure Owe consider as before T = 100 and we plot the evolution of the error in Energy 
conservation law along the discrete solution with h = 0.05 and r = 0.001. 

Test 6.2. We consider problem ([5|) with material tensor, load and initial data given by 

]K = / 2 , /(m) = sin(it), mo( 2;) = 0, uo(x) = sin(7ra;i) sin(7ra;2). 

We implement the fully discrete problem in the time interval [0,1] with the SIM coupled with 
the MFD discretization for the usual sequence of polygonal meshes introduced above. In Table 
?? we choose r = 0.001 and we show the errors for the Hamiltonian ah,T and for the Energy 
eh,T for different values of the mesh size h. We observe that we achieve the theoretical order 
h? of convergence. In Table ?? we fix the mesh size h = 0.05 and we display the errors for the 
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Figure 5: Behaviour of discrete Hamiltonian functional along the sequence with h = 

0.05 and r = 0.001. 



Figure 6: Behaviour of Energy conservation law error along the sequence {u^, with h = 0.05 
and r = 0.001. 

Hamiltonian and for the Energy as a function of t. In this case we observe that the error for 
the Hamiltonian is almost constant in r: the error due to the spatial discretization dominates 
the time component of the error. 



Table 2: Sh^r 

and eh,T for fixed time step size r = 

0.001. 


h = 0.2 

h = 0.1 

h = 0.05 

h = 0.025 

^h.T 

7.5726593e - 03 
2.2022797e-02 

1.9485268e-03 
7.5684283e - 03 

5.0100734e-04 

1.3466833e-03 

1.2486786e-04 
3.3563875e- 04 



Table 3: Sh,T 

and eh,T for fixed mesh diameter h 

= 0.05. 


T = 0.1 

T = 0.05 

r = 0.025 

T = 0.0125 

^h.T 

^h,T 

4.8131147e-04 

2.2022797e-02 

4.9593346e - 04 
7.5684283e - 03 

4.997307e - 04 
1.3466833e-03 

5.0068913e- 04 
3.3563875e - 04 


In Figure 0 we consider a larger final time T = 100 and we plot the behaviour of the error 
Sh,T in the discrete Hamiltonian functional along the sequence (u^, with respect to the initial 
value for h = 0.05 and t = 0.001. 

In Figure [5] we consider again T = 100 and we show the evolution of the error in Energy 
conservation law along the discrete solution with h = 0.05 and r = 0.001. 


7 Conclusions 

In this paper we have analysed the structure and the time invariants of the nonlinear wave 
equation discretized by mimetic approach. We have proved that the MFD discretization pre¬ 
serves the hamiltonian formulation of the problem and that the Hamiltonian and the Energy 
are still semi-discrete invariants of the solution. We have also derived a convergence theory for 
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Figure 7: Behaviour of discrete Hamiltonian functional along the sequence with h = 

0.05 and r = 0.001. 



Figure 8: Behaviour of Energy conservation law error along the sequence (rtJJ, v^) with h = 0.05 
and T = 0.001. 


the method, obtaining an order for the discrete norm of the error among the solution 
of the continuous and discrete problems. We have then considered the fully discrete scheme 
by making use of the MFD method coupled with the SIM time integrator: we have derived 
the convergence rate of the method and we have investigate the behaviour the Hamitonian and 
Energy. In light of these results we belive that the spatial discretization by making use of the 
MFD technique is a good choice in the context PDFs with conservation laws. In the present 
manuscript the focus is on the spatial discretization, thus the use of the mid-point scheme, for 
the time discretization, should be simply understood as a model symplectic method. 
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